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The identification of coherent structures from experimental or numerical data is an 
essential task when conducting research in fluid dynamics. This typically involves the 
construction of an empirical mode base that appropriately captures the dominant flow 
structures. The most prominent candidates are the energy-ranked proper orthogonal 
decomposition (POD) and the frequency ranked Fourier decomposition and dynamic 
mode decomposition (DMD). However, these methods fail when the relevant coherent 
structures occur at low energies or at multiple frequencies, which is often the case. 
To overcome the deficit of these “rigid” approaches, we propose a new method termed 
Spectral Proper Orthogonal Decomposition (SPOD). It is based on classical POD and it 
can be applied to spatially and temporally resolved data. The new method involves an 
additional temporal constraint that enables a clear separation of phenomena that occur 
at multiple frequencies and energies. SPOD allows for a continuous shifting from the ener¬ 
getically optimal POD to the spectrally pure Fourier decomposition by changing a single 
parameter. In this article, SPOD is motivated from phenomenological considerations of 
the POD autocorrelation matrix and justified from dynamical system theory. The new 
method is further applied to three sets of PIV measurements of flows from very different 
engineering problems. We consider the flow of a swirl-stabilized combustor, the wake 
of an airfoil with a Gurney flap, and the flow field of the sweeping jet behind a fluidic 
oscillator. For these examples, the commonly used methods fail to assign the relevant 
coherent structures to single modes. The SPOD, however, achieves a proper separation 
of spatially and temporally coherent structures, which are either hidden in stochastic 
turbulent fluctuations or spread over a wide frequency range. The SPOD requires only 
one additional parameter, which can be estimated from the basic time scales of the flow. 
In spite of all these benefits, the algorithmic complexity and computational cost of the 
SPOD are only marginally greater than those of the snapshot POD. 


1. Introduction and motivation 

1.1. Contemporary methods for data reduction 

Today’s high fidelity computational fluid dynamics (CFD) and high-end experimental 
data acquisition systems tend to produce vast amounts of data that are getting harder 
to interpret and overview. Methods to analyze such data are numerous and are always 
developing to stay in line with acquisition and computation systems. The most challeng¬ 
ing data stem from turbulent flows that feature a huge range of temporal and spatial 
scales. A key challenge in turbulent flow data mining is the distinction of deterministic 
coherent motion from purely stochastic motion. Numerous methods exist that exploit 
the periodicity or energetic dominance of these coherent structures. These methods 
range from classic Fourier decomposition to dynamic mode decomposition (DMD) and 
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proper orthogonal decomposition (POD). The most prominent among them are shortly 
introduced in the following. 

POD has been widely used since its introduction by Lumley (1970) and Sirovich (1987). 
It was applied in nearly every fluid dynamic field. Beyond fluid dynamics, this method is 
also known as singular value decomposition, principal component analysis or Kahunen- 
Loeve expansion (Berkooz et al. 1993). The basic idea behind this method is to construct 
an optimal basis that represents most of the data variance with as few basis functions 
as possible. In context of POD the variance is turbulent kinetic energy. Therefore, the 
POD searches for the most energetic modes whereby coherent structures with high energy 
content are likely to be represented by POD basis functions (Holmes et al. 2012). 

Another classical approach is the linear stochastic estimation introduced by Adrian & 
Moin (1988), where the readings of different sensors are related via a linear mapping. This 
is closely related to the extended POD (Boree 2003), also described in a unified framework 
(observable inferred decomposition) by Schlegel et al. (2012). In recent extensions of linear 
stochastic estimation, the use of time delays between the different sensors and also the 
use of one sensor at multiple time instances is pursued to separate periodic coherent 
structures from turbulent fluctuations (Durgesh & Naughton 2010; Lasagna et al. 2013). 
This approach was also used to improve the determination of harmonic POD modes 
from few pressure sensors (Hosseini et al. 2015). These utilisations of data from various 
time instances are also related to the temporal constraint used for the POD extension 
proposed in this article. 

Targeting the temporal periodicity of the coherent structures, spectral methods like 
discrete Fourier transform (DFT) and the recently introduced DMD (Rowley et al. 2009; 
Schmid 2010) come into play. These methods commonly span the mode space according 
to fixed frequencies, which enables the identification of coherent structures within small 
spectral bandwidths. In contrast to the DFT, the DMD also distinguishes modes with 
respect to their linear amplification. The recently introduced extended DMD (Williams 
et al. Submitted) tries to overcome the limitations encountered by the (linear) DMD 
approach when trying to decompose data from nonlinear systems. The idea is to use 
nonlinear functions that create observables of the data, which are exactly described by 
a linear system. This approach translocates the problem towards the identification of 
these nonlinear functions, which can be solved using the “kernel trick” that is common 
in machine learning. This paper presents an alternative approach, which extends POD 
to account for temporal dynamics in addition to energetic optimality. 

1.2. Why yet another method? 

After this short and incomplete review of data processing methods, one may ask if 
there is need for another method. The answer is probably no, so we take the most used 
method (POD) and bring it up to date for present research issues. The approach pursued 
here includes a simple yet effective extension to the classical POD, which leads to a 
more general method comprising POD and also the DFT. This approach unifies existing 
methods, but also offers possibilities beyond these. From the authors’ experience, the 
currently available methods often fail when applied to challenging flow data. These stem 
from flows with weak coherent structures where the recorded data have low signal to 
noise ratios, from flows with intermittent dynamics, or from flows featuring multi-modal 
interactions leading to frequency modulations, to name a few. In such cases, much effort 
is required to optimize the data processing until satisfactory results are obtained. The 
usual escape route is to focus on a certain spatial region or to apply suitable filters 
to pick out a certain wavelength or frequency range. This involves trial and error or 
requires prior knowledge of the investigated flow. There is also the danger of cutting 
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off a substantial portion of the data, leading to false interpretations. These procedures 
can be collected under the heading “identifying symmetries” as done by Holmes et al. 
(2012). The drawback of this approach is that the investigated flow must feature some 
symmetries and they must be known a priori. A recent application shows the huge variety 
of spatial and temporal filtering together with POD to separate different phenomena into 
different modes (Bourgeois et al. 2013), exemplifying the complexity of this approach. 

The usage of spectral methods for highly turbulent flows is even more challenging than 
POD. The variable frequency of single coherent structures and intermittent occurrence of 
different structures with the same frequency hinders a proper decomposition. In terms of 
the DPT, averaging of spectra from multiple measurements or sensors is essential to get 
reliable results. Analogously, for the DMD, averaging over several events is an option to 
reject noise (Tu et al. 2014). Nonetheless, the results obtained with DFT and DMD suffer 
from limiting the temporal dynamics to single frequencies. Turbulent flows hardly ever 
feature discrete frequencies and it is not always valuable to restrict a single mode (flow 
phenomenon) to a single frequency. Coherent structures that feature signihcant phase 
jitter or frequency modulation are represented by many modes at similar frequencies. In 
contrast, the POD puts no temporal constraint on the modes. This can result in modes 
that represent flow phenomena occurring at largely different temporal scales. Thus, it is 
often hard to interpret these modes and draw meaningful conclusions from the temporal 
dynamics. 

From our point of view, there is a big gap between the energetically optimal decompo¬ 
sition of POD and the spectrally clean decomposition of DFT or DMD. This gap will be 
bridged with the spectral proper orthogonal decomposition (SPOD) introduced in this 
article. This new method not only places itself somewhere in between these two extrema, 
but it allows for a continuous shifting from one to the other. The main idea is to apply a 
filter operation to the POD correlation matrix, which will force the POD towards clear 
temporal dynamic. Depending on the filter strength we continuously sweep from classic 
POD to DFT. 

The remainder of this article is organized as follows: The proposed method is described 
in detail in §2. The reader is guided from snapshot POD via an in-depth interpretation 
of the correlation matrix towards the general description of the SPOD. In addition, a 
method is explained to identify coupled mode pairs describing a single coherent structure, 
which becomes handy when working with SPOD. In §3, the new method is demonstrated 
on three different experimental data sets. The results are compared against POD and 
DFT to point out the benefits of SPOD. In §4 the capabilities of SPOD are summarized, 
based on the findings from the application to experimental data. 


2. Description and interpretation of the proposed method 

2.1. Classical snapshot POD 

To introduce the method and the nomenclature, the snapshot POD approach is 
described hrst. We start off with a decomposition of a data set into spatial modes and 
temporal coefficients: 

N 

U{x, t) = u{x) + u{x, t) = u{x) + ^ ai{t)(I>i{x) . (2.1) 

i=l 

Note that only the fluctuating part u{x, t) is decomposed. It is split into a sum of spatial 
modes and mode coefficients Oi. A set of M spatial points recorded simultaneously over 
N time steps is considered. To calculate the POD, the correlation matrix of this data set 
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is needed. For data obtained from particle image velocimetry (PIV) or CFD, the number 
of spatial points is usually larger than the number of snapshots. The correlation matrix 
is then calculated between individual snapshots (temporal correlation). The alternative 
approach (spatial correlation) that applies for M ^ N, is detailed in appendix A. The 
correlation between two snapshots is calculated from an appropriate inner product (,), 
usually dehned as the inner product 

{u{x),v{x)) = j u{x)v{x)dV, (2.2) 

y 

where V specifies the spatial region or volume over which the correlation is integrated. 
The elements of the correlation matrix R are given by 

^ {u{x, t,),u{x, tj)). (2.3) 

Matrix R is of size N x N. 

The temporal coefficients ai and mode energies Ai are obtained from the eigenvectors 
and eigenvalues of the correlation matrix. 

RcLj^ = Xidi i Al ^ A 2 ^ * * * ^ Aiv ^ 0 (^•^) 

The subscript i refers to single eigenvalues, which are sorted in descending order. Since 
the di are the eigenvectors of the real symmetric positive-definite matrix R, they are 
orthogonal. Moreover, they are scaled with the energy of the single modes such that 

-^(cti,Q-j) — XiSij^ (^•^) 

where (,) denotes the scalar product. The spatial modes are obtained from the projection 
of the snapshots onto the temporal coefficients 

1 ^ 

Mx) = ( 2 . 6 ) 

* 

These modes are orthonormal by construction 

(<l>.,<l>,)=%. (2.7) 

The formulation so far is perfectly in line with classical snapshot POD, which can also 
be computed by a singular value decomposition. However, since the SPOD requires a 
manipulation of the correlation matrix we retain the classical form. 

2.2. Properties of the correlation matrix 

The SPOD described in this article is essentially a filter applied to the correlation ma¬ 
trix R. To offer a better understanding of this approach, the structure of the correlation 
matrix R is inspected Hrst. 

Figure 1(a) shows the structure of the correlation matrix for the data set of a forced 
turbulent jet. The data were acquired with PIV inside a 2D-plane aligned with the jet 
axis. The considered flow shows strong vortex shedding at the forcing frequency (the 
acquisition frequency is 25 times the forcing frequency). The presence of these periodic 
patterns in the flow, and their convection within the observed flow field, lead to a diagonal, 
wave-like structure of the matrix. This is closely related to the periodicity of the auto¬ 
correlation coefficient. In fact, if the individual elements of the correlation matrix R are 
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Figure 1. (a) Pseudo-color plot of the correlation matrix elements and (b) the 

corresponding correlation coefficient R. The displayed data are picked from PIV measurements 
of a forced turbulent jet. 


summed up along the diagonals, we get the spatially averaged auto-correlation coefficient 


R{t) 


JJ u(x,t)u(x,t + T)dxdt 
JJ u'^{x,t)dxdt 


( 2 . 8 ) 


It is depicted in figure 1(b), showing the same periodicity as the correlation matrix. The 
auto-correlation coefficient itself represents the spectral content of different timescales 
and wavelengths and it is directly related to the power spectral density of the underlying 
data. However, it contains no information of the phase of individual frequencies, due to 
the reference of the signal to itself. This is why the elements along the diagonals of R 
look so similar, as they represent only relative changes with respect to the time step on 
the main diagonal. Thus, increased similarity along the diagonals of R is equivalent to 
an increased similarity of the dynamics of the underlying signal. This property will be 
discussed more deeply in section 2.4. The obvious consequence from these hndings is: 
If we want to obtain smooth dynamics from the POD, we have to enforce the diagonal 
similarity of the correlation matrix. This is where we step into spectral POD. 


2.3. General description of the SPOD 

The yet so simple, but radical approach is a filter operating on the correlation matrix 
R. To augment the diagonal similarity of R a simple low-pass filter is applied along the 
diagonals. This results in a filtered correlation matrix S with the elements given as 

Nf 

~ ^ ^ PkPi+kJ+k- (^* 9 ) 

k=-Nf 

The filter above is just a symmetric finite impulse response filter with a filter coefficients 
vector g of length 2Nf -b 1. The most simple approach would be a box filter, where 
all coefficients have the same value gk = 2 n]+i ■ examples discussed later, we 

use a Gaussian filter, which features a smooth response in time and frequency domain. 
Moreover, we choose a standard deviation as such that the filter gives the same cut-off 
frequency as a box filter with half the length. In fact, any kind of digital finite impulse 
response filter can be used here. 

The further procedure of the SPOD is the same as for the classical POD. From the 
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correlation Matrix S the temporal coefficients bi and mode energies fii are obtained from 
the eigen decomposition 

Sbi=^ibi\ Ml ^ > • • • ^ ^ 0. (2.10) 


The temporal coefficients are also scaled with the mode energy and they are still 
orthogonal 

— {bi,bj)=iii5ij. ( 2 - 11 ) 

The spatial modes are hnally obtained from the projection of the snapshots onto the 
temporal coefficients 


N 


^^{x) = 


N^i, 


^ ' bi{tj)u{x, tj), 


( 2 . 12 ) 


where these modes are no longer orthogonal. This property of the spatial modes is detailed 
in appendix B. The total energy of the data set is still represented by the decomposition 
Xi = X) M*)) 111® energy per mode is less for the first modes. Hence, increasingly 

plain temporal dynamics are obtained at the expense of spatial orthogonality and a 
dispersed SPOD spectrum. Nevertheless, the decomposition (as in (2.1)) 


N 

U{x,t) =u{x)+'^bi{t)'l't{x), (2.13) 

i=l 

is still exact if all N SPOD modes are used for the re-composition. 

If the filter size is extended over the entire time-series, the filtered correlation matrix 
converges to a symmetric Toeplitz matrix. This matrix has the form 

S,,,=R{At\i-j\), (2.14) 

with the diagonals given by the average correlation coefficient (2.8). This special matrix 
is also known as the covariance matrix and its eigenvalues are tracing out the power 
spectral density of the underlying time series (Wise 1955). This equality is a part of Szegos 
theorem and it is valid for the limiting case where the number of samples approaches 
infinity. To discuss this feature for finite series, the treatment of the start and end of 
the time series must be clarified. At the boundaries of /?, the filter operation is not 
properly defined, since the symmetric filter lacks elements before and after the finite 
series. These elements can either be replaced by zeros or the time-series is assumed to 
be periodic. In the latter case, a box hlter of the same size as the number of snapshots 
delivers a symmetric circulant matrix, with its eigenvalues and eigenvectors given by the 
Fourier transform of the first row (Gray 2005). The DFT and the SPOD produce the 
same decomposition for this limiting case. Hence, the SPOD is able to continuously fade 
from the energetically optimal POD to a purely spectral DFT. What happens in between 
these two limits is very promising, as will be shown later. 


2.4. The correlation matrix from a dynamic systems point of view 

To further consolidate the previous considerations, the correlation matrix is represented 
in an analytical framework. For the moment, the temporal evolution of the investigated 
flow is assumed to be locally governed by a linear time invariant model. The term locally 
refers to a short and finite temporal extension, and it is closely related to the hlter size 
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Nf. The temporal evolution of this system is given by 

^=Lu{t), (2.15) 

where the matrix L is the system matrix describing the entire dynamics and the spatial 
points of the velocity filed are organized as rows in u. Starting from a reference snapshot 
Uq = u(t = 0), the field at any time step can be calculated by the matrix exponential 

u(t) = e^*Uo. (2-16) 


To allow further simplifications, we require the matrix L to be normal. This allows for 
the decomposition L = UDU*, where t/ is a unitary matrix (eigenvectors of L), D is 
a diagonal matrix (eigenvalues of L) and * means the conjugate transpose of a matrix 
(adjoined). Hereafter, equation (2.16) becomes 

u{t) = e^‘^^'*Uo = Ue°^U*Uo. (2.17) 

The diagonal elements of the matrix D contain the complex eigenvalues dj- of the system. 
Each of these eigenvalues contain the amplification rate a and the frequency w of the 
related mode dk — (Jk + ioJk (i = V—1). The diagonal matrix can be decomposed in 
amplihcation rate L and frequency O, thus D = L + iO. Note that for this linear 
approach, the non-linearity and non-normality of the Navier-Stokes equations contribute 
to parameter variations of the linear system. The flow is assumed to behave like a 
linear normal system within the time scale of the filter and all non-linear and non¬ 
normal dynamics are represented by variations in a and w. This kind of approach is 
also pursued in the generalized mean field model of Luchtenburg et al. (2009), where the 
mode interaction via the mean flow is represented by an interaction of linear oscillators 
trough nonlinear coupling of model parameters. 

With the formulation in (2.17), the inner product, which forms the elements of the 
correlation matrix, is simplified to 


'u{x,ti),u{x,t 2 )) = Ue^^^U*Uo) 

(2.18) 

= (^U*Uo, (ye°‘i)* Ue^*W*uo) 

(2.19) 


(2.20) 


(2.21) 


(2.22) 


An inspection of the exponent in (2.22) reveals that the inner product only depends on 
the sum and difference of time steps. According to this equation, the velocity fields are 
projected onto the subspace spanned by the linear operator U*Uo and changes of the 
inner product are only governed by the eigenvalues of the system. 

Within the context of the correlation matrix, we use the abbreviated nomenclature 
for the snapshots Ui = u{iAt) and for the projected velocity u = U*Uq. The correlation 
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matrix constructed around the neighborhood of snapshots Uq, yields 



f (m_i,m_i) 

(mo,M-i) 

(mi,m_i) \ 


Ssub — j 

(m_1,Mo) 

(mo,mo) 

(mi,mo) 

(2.23) 



(mo,Mi) 

(Mi,Mi) / 


i 

(m, 

u) {u, 

(ft, \ 



= I {u,e (it, m) (it, ] . (2.24) 

y (m, e^'^^*it) (m, (it, e^^^*ii) J 

The actual properties of this complex expression are highlighted by showing just the 
factors for L and iO in the exponents, given as 


1 

to 

-1 

0 ■ 


■ 0 

-1 

1 

to 

-I 

0 

I 

At; iO: 

I 

0 

1—1 

1 

1- 

o 

1 

to 

1_ 


_1 

1 

-1 

o 


(2.25) 


It is perfectly visible that the changes of the correlation matrix along the diagonals are 
caused by amplification of modes, while changes along the anti-diagonals are caused 
by the modes’ frequency. Therefore, the changes of the correlation matrix along the 
diagonal are related to the modes amplitude and the crosswise changes are caused by the 
modes phase. This was already observed in section 2.2, where the correlation matrix was 
compared to the temporal auto-correlation. Note that the diagonals in the figures start 
at the lower left corner whereas in the formulas they start in the upper left. 

As mentioned in section 2.3, the SPOD and the Fourier transform become equal when 
the filter operation is extended over the entire sequence and periodic boundary conditions 
are applied. The correlation matrix becomes a Toeplitz (circulant) matrix for this case, 
and also the dynamics of the entire sequence are described by a linear system. The latter 
point is evidenced by equation (2.22). For the entire time series to be governed by a 
linear system, the complete correlation matrix must have the form of the subset shown 
in equation (2.24). Since there can’t be any amplification in case of periodic boundary 
conditions, the correlation matrix is solely defined by the frequencies of the modes (O). 
Thus, the matrix is constant along the diagonals, which is the presumed shape. 

Recall that the SPOD filter operation acts along the diagonals of the correlation matrix. 
The excursion to system dynamics indicates that the SPOD puts a constraint on the 
temporal variation of the mode coefficients. To clarify these properties we describe the 
coefficient and its derivative by an analytical signal 


b{t) = A(t)e‘^(‘) ; = [a{t) + icc(t)] b{t). (2.26) 

In this framework, changes along the diagonal of the subset (2.25) are caused by temporal 
variations of A or to. This behavior is sketched in figure 2, where the corresponding 
changes of the correlation matrix anti-diagonal are depicted. There, changes of the 
amplitude cause an overall scaling and changes of the frequency move the zero crossings. 
The SPOD filter (2.9) smooths along the diagonals and therefore it equalizes the anti¬ 
diagonals of consecutive time steps. Thus, the two curves in figure 2 move closer together, 
which limits the rate of change of the amplitude dAjdt and the frequency dcojdt. 

The exact spectral properties of the mode coefficients might not be directly concluded 
from the considerations in this section, bnt they become obvious in the SPOD results. 
They can be summarized as follows: The low-pass filter applied in equation (2.9) defines 
a certain spectral bandwidth. The filter response and cut-off frequency can be obtained 
from the coefficients gi. The individual SPOD coefficients bi feature low-pass filtered 
amplification rates cr and band-pass filtered frequencies w. For higher SPOD modes the 
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correlation matrix 


Ri, 



change of A 


\\ change of w 


. tk 

- tk S' At 


I* -il 


Figure 2. Schematic that illustrates the influence of parameter changes on the correlation 
matrix anti-diagonals. The correlation matrix on the left indicates the orientation of the 
diagonals (different times t) and the anti-diagonals (different time delays t = \i — j\At). The 
parameters A and uj are specified in (2.26). 


effect of the filter becomes less pronounced. This is plausible, since the filter is applied 
to the entire correlation matrix and not as a constraint for the single modes. 

2.5. Identification of coupled modes 

One crucial point in POD and SPOD is the identification of linked modes. Assuming 
the presence of periodic coherent structures, their dynamics are described by a pair of 
modes, analogue to the sine and cosine in the Fourier space or the real and imaginary part 
of DMD modes. They constitute the real and imaginary part of the analytical coefficient 
in equation (2.26). The coupling of such a mode pair is not given by the SPOD and 
it has to be identified a posteriori. Coupled modes typically show a similar amount of 
energy and pair in the POD spectrum (Oberleithner et al. 2011, 2014). For more complex 
dynamics with multiple energetic modes, the pairs are not easily identified and visual 
inspection of Lissayous figures and spatial modes is required. This manual procedure is 
cumbersome and by no means objective. 

To provide an alternative, we propose an unbiased approach that gives a quantitative 
measure of the dynamic coupling of individual modes. Based on a DMD of the temporal 
coefficients, we derive a relation of the spectral proximity of the SPOD modes. The 
general procedure is schematically outlined in figure 3. From the depicted operations, the 
DMD and the mode coupling are discussed in this section. The Fourier decomposition of 
the coefficients may similarly provide a basis to compare the coefficients, but the DMD 
turned out to be more reliable for this task. 

For the DMD of the SPOD coefficients, their temporal evolution is assumed to be 
governed by a linear operator T 

b{t +At) = Tb{t). (2.27) 

To approximate this operator, the SPOD coefficients are arranged in two matrices X = 
[b(0) b{At) ... b{{N — 2)At)] and Y = [b{At) b{2At) ... b{{N — l)At)] (following the 
notation of Tu et at (2014)). Hereafter, the operator is given by 

T=YX-\ (2.28) 

where X^^ is the Moore-Penrose pseudoinverse of X. Alternatively, this can be solved 
as a least squares problem, minimizing ||TX— K||. To reject measurement noise in the 
identification procedure, only the “physical part” of the SPOD modes is considered for 
the calculation of the operator T. That means, only the modes with acceptable signal to 
noise ratio should be considered. Therefore, the number of retained modes is calculated 
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time-resolved 
snapshot sequence 


u{x,y,t) 





Figure 3. Schematic illustrating the main steps towards the identification of coupled modes 
(red and blue lines indicate real and imaginary parts of an analytic signal). The data displayed 
here were derived from measurements of a forced turbulent jet (see also figure 1). 


from the energy resolved by the SPOD, truncated after Nc modes, with 

£{Nc) = (2.29) 

In the examples shown later the modes are truncated around £{Nc) = 0.95. This value 
depends on the signal to noise ratio of the considered measurement, which can be 
estimated from the POD spectrum (Raiola et al. 2015). Note that the number of retained 
modes increases for wider SPOD filters and corresponding flatter SPOD spectra (/i_,). 
The DMD modes are obtained by the eigen-decomposition of matrix T as 

Tci = iy^c^. (2.30) 

The eigenvalues comprise the frequencies uJi and amplification-rates (7^ of the operator 
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T and are given by the logarithm of the eigenvalue ln(r'i)/At = CTi + iwi. The eigenvectors 
Ci hold the relative spectral content of the single SPOD coefficients with respect to these 
frequencies. More precisely, the element Cij of vector q is the spectral content of the 
single mode coefficient bj with respect to Vi. The actual modal representation is given by 

(2.31) 

i=l 


It must be noted that this decomposition is only exact if Nc = N, whereas in the current 
approach Nc < N. Nevertheless, the decomposition gives a common spectral basis, which 
allows the ranking of spectral similarity of the temporal coefficients b{t). The developed 
proximity measure is given by 


Cij = Im 


Nr 


'^Ck,icl jSgn (Im(i^fe)) 


.fc=i 


(2.32) 


where the coefficients are normalized to ( 0 ^, 0 ^) = 1. The sign (sgn) in this expression 
accounts for the conjugate pairs that appear in the DMD spectrum (mirrored at the real 
axis). 

For two modes to be coupled, they must have a similar spectral content, which is 
either shifted a quarter period forward or backward. This implies a positive or negative 
imaginary part of the harmonic correlation (2.32), respectively, and coupled modes appear 
as peaks in the matrix C. Hence, the row and column indices of the maximum of C identify 
the first coupled SPOD modes. The corresponding row and column in C are then set 
to zero and the next lower maximum is identified. This procedure is repeated until all 
modes are paired. It has to be noted that this approach also creates weakly correlated 
and possibly unphysical mode pairs. 

Together with the identification of coupled modes, the procedure gives an average fre¬ 
quency of the coherent structure represented by the mode pair. Therefore, the eigenvalues 
i^k of the matrix T are sorted in descending order with respect to their content for the 
identified mode pair Ck = ^ j. The frequency is given as the weighted sum of the 

eigenvalues 


ELiIm{ln(:y,)}gi 
2'KAt Yh=i D 


(2.33) 


The weighting accounts for the relative energy content of a mode pair with respect to 
the single frequencies. In fact, only the most relevant eigenvalue (n = 1) can be picked 
to determine the frequency, but for practical application it is recommended to use more 
than one eigenvalue as noise may corrupt them. For the examples discussed in the next 
chapter, we used an average over three eigenvalues (n = 3) to get accurate results. 

The coupled SPOD modes are considered as one complex mode (see equation (2.26) 
and figure 3) similar to the Fourier mode. The relative energy content of the identified 
modes is computed as 


+ hj 

~ ^Nr ’ 


(2.34) 


where i and j again refer to the indices of the coupled SPOD modes. 
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3. Applications to experimental data 

In this chapter the SPOD is applied to three different data sets. All three examples 
originate from very different engineering problems, demonstrating the capability and 
broad applicability of SPOD. We consider the flow of a swirl-stabilized combustor, the 
wake of an airfoil with Gurney flap, and the flow field of a sweeping jet generated with a 
fluidic oscillator. All three flows were recorded with the same PIV measurement system. 
It consists of a Photron Fastcam SA 1.1 high-speed camera (IMpixel at 2.7kHz double 
frame) and a Quantronix Darwin Duo laser (30mJ at IkHz). The PIV data were processed 
with PIVview (PIV TEC GmbH) using standard digital PIV processing (Willert & Gharib 
1991) enhanced by iterative multigrid interrogation with image deformation (Scarano 
2002),(Raffel et al. 2007, pp. 146-158). Analyzing the present data sets with existing 
POD, DFT or DMD approaches caused some difficulties. As will be demonstrated, the 
SPOD is able to handle these shortcomings. The DMD and the DFT equally suffer from 
the restriction of the modes to narrow frequency bands, therefore we limit the following 
presentation to DFT. This choice is particularly handy as the DFT is a limiting case of 
the SPOD. 


3.1. Swirling jet undergoing vortex breakdown 

At Hrst, we consider the flow field of a swirl-stabilized combustor. Swirling jets are 
widely used in the gas turbine industry due to their capability of obstacle-free flame 
anchoring and enhanced mixing. The experimental setup to study these flows is sketched 
in figure 4. Swirl is generated by injecting fluid tangentially into a mixing tube that 
terminates in the combustion chamber. The cylindrical-shaped chamber is made of quartz 
glass to allow optical access for PIV. Flow measurements are conducted in the meridional 
section as indicated in the schematic. The case investigated here is non-reacting at a 
Reynolds number of 58 000 based on the nozzle diameter D and the bulk velocity at 
the nozzle exit. Additional details about the experimental setup can be found in Reichel 
et al. (2015). 

The mean flow field and the spatially-averaged power spectral density (PSD) is depicted 
in figure 5, with the Strouhal number based on the same length and velocity scale as the 
Reynolds number [St = /D/rtbuik)- The flow exhibits a strong recirculation zone in the 
center, surrounded by an annular, strongly diverging jet. The cross-sectional jump at the 
combustion entrance creates an additional external recirculation zone between the jet 
and the confining walls. The spectral content of the flow is spread over all time scales 
and it decreases with increasing Strouhal number. The spectrum gives no indication of 
any dominant coherent structure even though these flows typically feature helical global 
modes (Oberleithner et al. 2011). 

Figure 6 (a), (b), and (c) illustrate the results from the SPOD for the filter lengths Nf = 
0, 10, and 2000, respectively. Note that the limiting cases Nf = 0 and Nf = 2000 produce 
results equivalent to those obtained with classical POD and DFT respectively, while 
the case in between represents the SPOD. Hence, this particular presentation concisely 
demonstrates the difference between POD, SPOD and DFT. 

Each of the three cases in figure 6 show the so-called SPOD spectrum, where every 
mode pair is represented by a single dot. The size and color of the dots indicate the 
harmonic correlation of a mode pair Cij, according to equation (2.32). The frequency of 
a mode pair is determined according to (2.33) and the energy from the two eigenvalues 
relative to the total energy from (2.34). On the right side of every case in figure 6, three 
spatial modes Ti(x) and the corresponding temporal coefficients bi{t) are posed above 
each other. The spatial modes are visualized by the crosswise velocity component (in 
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Figure 4. Experimental setup of the swirl stabilized combustor. Air enters from the left, passes 
a swirl generator and exits into the combustion chamber. Flow field measurements with PIV are 
conducted in the meridional plane as indicated by dashed square (ROI). 
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Figure 5. Swirl-stabilized combustor flow: Time-averaged flow field depicted by (a) contours 
of velocity magnitude and streamlines, and (b) spatially-averaged power spectral density. 


y-direction) together with streamlines of the time-averaged flow. They are numbered 
likewise in the SPOD spectrum and between the small mode plots. The plots are given 
without axis labels to allow a compact representation of the data, the section is the 
same as for the mean flow shown in figure 5(b). The time coefficients are represented by 
their power spectral density, where the time series is split into five (50% overlapping) 
sections, which are Fourier transformed and averaged. The horizontal dotted lines in the 
PSD plots indicate a spacing of three orders of magnitude (1000) and the i/-axis is scaled 
logarithmically. The spectral averaging was also applied for the power spectra shown in 
figure 5(b) and in the subsequent PSD plots. 

The POD (figure 6(a)) yields a broad spectrum of modes, where modes at lower 
Strouhal numbers have more energy. There are several modes with high harmonic 
correlation (diameter and color of the points), and high energy contents K. The spatial 
shape of the low-frequency mode (label 1) shows clear spatial symmetry and a limited 
spectral bandwidth {St « 0.1). This mode is frequently observed in swirl-stabilized 
combustors and it is associated with a global hydrodynamic instability (Terhaar et al. 
2015). From the four additional outstanding modes between St = 0.3 and St = 0.8 
we pick two for further investigation. Their spatial structures show no clear symmetries 
and indicate mixtures of several spatial wavelengths. Accordingly, the mode spectra are 
broad and show only a slight hump at the frequencies indicated by the SPOD spectrum. 





























































Figure 6. Swirling jet: Results from SPOD for different filter lengths (a) Nf = 0 (POD), (b) 
Nf = 10 (SPOD), and (c) Nf — 2000 (DFT). For every filter length the SPOD spectrum is 
displayed as scatter plot (left), where a single dot indicates one mode pair (size and color Cij 
in (2.32)). For three selected pairs the spatial modes (upper row) and PSD of the temporal 
coefficient (lower row) are depicted. They are indicated by numbers in the SPOD spectrum, as 
well as between the small mode plots. 






















































































































































Spectral proper orthogonal decomposition 15 

The other modes around St = 0.5, which are not shown here, show similar spatial and 
spectral content. 

Overall, the POD indicates the presence of a single mode at low frequency, together 
with other coherent structures that are not properly resolved. The first most energetic 
(not inspected here) corresponds to a low frequency, large wavelength fluctuation, as 
indicated by the SPOD. Such slow changes of the mean flow are usually named shift 
modes (Luchtenburg et al. 2009; Hosseini et al. 2015). In this particular case the shift 
mode stems from weak movements of the vortex breakdown bubble. 

Consider now the SPOD in figure 6(b), with a filter length Nf = 10; From the SPOD 
spectra we identify three peaks at St = 0.09, 0.5 and 0.8. The first mode is the same as 
the one already identified by the POD, but its spectral content at higher frequencies is 
reduced. It describes a single-helical structure in the wake of the recirculation zone. The 
second identified mode exhibits a broad spectral peak at St = 0.5. Its spatial structure 
and Strouhal number match the global mode identified by Oberleithner et al. (2011). It 
is a single-helical mode linked to the precessing motion of the recirculation zone. The 
spatial structure of mode three suggests a double-helical shape. It is not a harmonic of 
the second mode, as their frequencies are not related. 

When the filter size is extended to its maximum, we get the the DFT (figure 6(c)) 
and the SPOD spectrum converges to the averaged PSD (figure 5(b)). The temporal 
coefficients converge to sines and cosines and all mode pairs show full harmonic corre¬ 
lation (uniform dot size in the SPOD spectrum). Since the selection based on harmonic 
correlation is impossible, we resort to the frequencies already known from the SPOD 
at Nf = 10. The spatial structures resemble the ones from figure 6(b), but they are 
corrupted by noise. Moreover, the spatial symmetries are no longer as clear as for the 
Nf = 10 case. Note that the corresponding spectral peaks are broadened due to the 
averaging procedure, which is applied here only for consistency. 

From this first example, we can point out some striking features of the SPOD. The 
SPOD is able to separate coherent fluctuation from stochastic turbulent fluctuations even 
though they both have the same energy contents (see the SPOD spectrum in figure 6(c)). 
The classical POD yields partially mixed structures that cannot be assigned to distinct 
flow phenomena, whereas the SPOD properly separates these structures and identifies 
them from harmonic correlations. The DFT instead shows the same structures at the 
identified frequencies, but they are corrupted with noise and the method itself would 
give no clue about the frequencies of interest. 

The structures identified with the SPOD may also be found with the POD if the 
decomposition is applied to a subsection of the measured domain. Moreover, the ex¬ 
ploitation of spatial symmetries prior to the POD decomposition usually provides good 
results for this type of flows (Terhaar 2015). Nevertheless, these alternative approaches 
would require prior knowledge of the shape or spatial extent of the structures, whereas 
the SPOD requires none of these. 

All modes identified by the SPOD show clear spatial symmetries and distinct spectral 
peaks. The frequency and shape of the first mode coincide well with previous experimental 
observations in swirl-stabilized combustors (Terhaar et al. 2015). The second mode is very 
similar to the one observed in unconfined swirling jets (Oberleithner et al. 2011). However, 
the presence of these different modes in a single flow configuration raises the question 
about their physical relevance. To deal with this issue, we conducted a linear stability 
analysis of the underlying mean flow, following the procedure outlined by Oberleithner 
et al. (2015). This analysis similarly delivered three unstable modes whose frequencies 
and azimuthal and axial wavenumbers match the SPOD modes surprisingly well. To limit 
the scope of this paper the analysis is not further detailed here. One important parameter 
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of the SPOD is the filter size Nf, which is twice the period of the second mode. The 
experiences gained throughout the first application show that a filter size of one to two 
periods of the mode of interest gives the best results. 

3.2. Airfoil with Gurney flap 

The second flow configuration considered here is the flow behind a pitched airfoil 
equipped with a Gurney flap. The experimental setup is shown in figure 7. It illustrates 
the working principle of the Gurney flap deployed at the pressure side of the airfoil. The 
flap creates additional lift (and drag), which can be used to locally control varying loads 
on large wind turbine airfoils (Bach et al. 2015a, 2014). The flow features around the 
Gurney flap are characterized by a single vortex that develops upstream of the flap and 
periodic shedding in its wake. The vortex upstream of the flap continuously grows up 
to a critical size, then it is shed into the wake, and it starts growing again. Here, a FX 
63-137 airfoil at 5° angle of attack is investigated at a Reynolds number of 180 000 based 
on chord length. The reference length for the Strouhal number is the flap height, which is 
3.6% of the chord. The measured region comprises only the wake of the airfoil capturing 
the dominant vortex shedding. More details about the experimental setup can be found 
in a preceding publication of these data (Bach et al. 2015&). The Strouhal number in the 
following results is calculated with the flap height h and the free stream velocity. 

The mean flow shown in figure 8(a) reveals a velocity deficit in the wake of the 
Gurney flap, which generates the vortex shedding. The PSD (figure 8(b)) indicates strong 
oscillations at St = 0.105 with a weak higher harmonic. In a previous investigation it was 
found that the vortex, which is shed from upstream of the flap causes an alteration of the 
periodic vortex shedding behind the flap (Troolin et al. 2006). Hot-wire measurements in 
the wake of the Gurney flap supported this assumption. The combination of a strong 
periodic flow pattern and the intermittent short-time events provides a formidable 
benchmark for the SPOD. 

The presentation of the decomposition with the different methods is organized in the 
same way as for the previous example. The classic POD decomposition is shown in figure 
9(a). The vortex shedding is represented by the most energetic POD mode with the 
highest harmonic correlation. The remaining modes show weak harmonic correlations 
and no distinct peak in the SPOD spectrum. The modes labeled 2 and 3 exhibit a 
broad spectral content with a spatial extent limited to the vicinity of the flap. There 
are additional modes with similarly compact spatial extent located further downstream. 
These compact modes describe the intermittent alteration of the vortex shedding during 
the passage of the single vortex that was generated upstream of the Gurney flap. An 
inspection of their time coefficients (not shown) reveals that these modes are only active 
one after another during one shedding period. Depending on the phase lag between the 
natural oscillation and the shedding of the upstream vortex, the developing wake vortex 
is either strengthened or weakened. The convection of this altered vortex is described by 
the spatial series of modes. This behavior is indicated by a complex interaction of these 
modes and the periodic shedding modes, which is hard to identify in the POD expansion. 

The SPOD yields a much clearer set of modes (figure 9(b)). In addition to the shedding 
mode, the SPOD also uncovers three other modes, which are offset from the rest. The two 
modes that appear at similar frequencies capture the alteration of the vortex shedding 
during the passage of the single vortex. Their mode shape is similar to the shedding mode, 
but with larger spatial wavelengths and lower frequencies (see mode 2 in figure 9(b)). 
The interaction of the upstream vortex with the vortex shedding increases the vortex 
size and thus the wavelength in the wake. Assuming a constant convection speed in the 
flow, this mode consequently has a lower frequency. In case of the SPOD the alteration 
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Figure 7. Schematic of the airfoil equipped with a Gurney flap at the trailing edge. Streamlines 
indicate the surrounding flow and the vortex upstream of the flap. The measured section (ROI) 
is a streamwise cut in the wake of the airfoil, capturing the periodic shedding behind the flap. 




Figure 8. Wake of the airfoil with Gurney flap: Time-averaged flow field depicted by (a) contours 
of velocity magnitude and streamlines, and (b) spatially-averaged power spectral density. The 
origin of the coordinate system is located at the trailing edge. 


of the vortex shedding is captured by a single mode (pair) and is thus much easier to 
interpret. The third mode represents the second harmonic of the vortex shedding with 
a clear spectral peak and clean spatial mode with twice the wavelength of the shedding 
mode. This higher harmonic is completely missed in the POD. The SPOD filter size Nf 
is equivalent to three shedding periods, which is approximately equal to the traveling 
time through the measurement domain. 

The DPT shown in figure 9(c) reproduces the spectrum shown in figure 8(b). The 
natural mode and its higher harmonic can be identified from the spectral peaks. The 
corresponding mode shapes are similar to the SPOD, although the higher harmonic is 
corrupted with noise, resulting in a fragmented spatial mode. The DPT at the frequency 
of the second SPOD mode gives no indication of the structure identified before and 
the vortex-vortex interaction is completely missed. This is attributed to the fact that 
this phenomenon is highly intermittent with varying frequencies and amplitudes, which 
cannot be represented by a single-frequency mode. The same dilemma applies for the 
DMD. 

Por this example, the SPOD has shown its ability to separate dynamics with similar 
spatial structures and frequencies but very different energies. The spectral proximity 
and spatial similarity of the involved dynamics impede the application of POD. The 
modulation of the natural vortex shedding is represented by a natural mode with several 
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Figure 9. Airfoil with Gurney flap: Results from SPOD for different filter lengths (a) Nj = 0 
(POD), (b) Nf = 15 (SPOD), and (c) Nf = 2000 (DFT). For every filter length the SPOD 
spectrum is displayed as scatter plot (left), where a single dot indicates one mode pair (size 
and color Ci,j in (2.32)). For three selected pairs the spatial modes (upper row) and PSD of 
the temporal coefficient (lower row) are depicted. They are indicated by numbers in the SPOD 
spectrum, as well as between the small mode plots. 
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intermittent modes. The DFT, however, with its single frequency modes does not capture 
the modulation of the shedding at all. The frequency constraint imposed by the SPOD 
is sharp enough to split the natural shedding from the modulation and soft enough to 
allow for frequency and amplitude variations. Hence the SPOD again gives easy access 
to dynamic features of the flow, which cannot be found with other common methods 
of similar algorithmic complexity. There may be feature tracking approaches capable of 
identifying the dynamics in this case, but they usually require more computational effort 
and might not be as versatile as the SPOD. 

3.3. Fluidic oscillator 

In this example, SPOD is applied to the flow field of the sweeping jet generated from 
a fluidic oscillator. This device is essentially a nozzle with feedback channels, which 
cause a self-sustained oscillation of the jet. Figure 10 shows the approximate geometry 
of this device and indicates the meandering shape of the sweeping jet. The shape and 
motion of the jet resemble a traveling wave. These devices are used for active flow control 
applications, where the sweeping motion of the jet allows a much wider actuator spacing 
resulting in less energy consumption (Woszidlo et al. 2014). The data presented here stem 
from an experimental setup investigating the spreading and entrainment of sweeping jets 
(Woszidlo et al. 2015; Ostermann et al. 20156). The data are recorded at a Reynolds 
number of 37 000 based on the nozzle diameter D and the mean velocity in the nozzle. 
These scales are also used for later calculation of the Strouhal number. The mean velocity 
in figure 11(a) show that the PIV domain is moved off the jet center towards the lower 
part of the jet. Data points closer than x/D = 2 were distorted due to laser light 
reflections. The spectral content averaged over the PIV domain (figure 11(b)) shows 
a narrow dominant peak and at least five higher harmonics. The narrow peaks indicate a 
stable operation at the fundamental frequency, while the additional peaks suggest more 
complex dynamics. The key challenge of this data set is to accurately reconstruct the 
sweeping jet dynamics from a truncated measurement domain. 

Figure 12 shows the results from the SPOD for filter lengths Nf = 0, 30, and 2000. 
As in the foregoing examples, these filter setting span the range between both limiting 
cases (from the POD to the DFT). The spectrum attained with the POD (figure 12(a)) 
reveals distinct modes at the fundamental frequency of the oscillator (labeled as 1) and 
at higher harmonics. The mode at the third harmonic frequency (labeled as 2) shows 
a surprisingly high harmonic correlation. The PSD of the mode coefficients reveal that 
each mode is not limited to a single frequency. The additional peaks in the PSD are 
partly attributed to the fact that only part of the jet is measured. During one oscillation 
period, the jet leaves and enters the measurement domain, creating sharp transitions in 
the time domain and thus a series of higher harmonics in the frequency domain. Due 
to the purely statistical POD approach, these higher harmonics appear in every mode 
coefficient, which contradicts the idea of a proper modal decomposition. The mode that 
seems to represent the fifth harmonic in the SPOD spectrum (labeled as 3) shows no 
distinct peak at all in the PSD of the coefficient. Thus, the POD of this data set does 
not provide a proper separation of the fundamental and higher harmonic contributions. 

If the SPOD is applied instead, the fundamental and harmonic modes line up perfectly 
(figure 12(b)). Now, the harmonics are separated clearly up to the seventh harmonic. The 
spectral content and spatial shape are further examined for the fundamental, the third 
and seventh harmonic. The PSDs of the mode coefficients reveal narrow spectral bands. 
The corresponding mode shapes show an appropriate spatial symmetry, although the 
PIV domain is cropped shortly above the symmetry line. It is worth mentioning that the 
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Figure 10. Schematic of the experimental setup with the fluidic oscillator. Air enters from 
the left, passes the oscillator and exits into the unconfined ambient air. The angle of the jet 
leaving the oscillator sweeps periodically up and down. The measured region (ROI) captures 
the meridional plane of the jet’s near field. The oscillator has a square nozzle, hence the thickness 
of the jet normal to the plane is also D. 




(b) 


Figure 11. Fluidic oscillator: Time-averaged flow field depicted by (a) contours of velocity 
magnitude and streamlines, and (b) spatially-averaged power spectral density. 


broad peak in the PSD of the seventh harmonic indicates considerable frequency jitter, 
while the mode shape remains remarkably smooth and symmetric. 

The results obtained with the DFT are presented in hgure 12(c). The peaks in the 
SPOD spectrum clearly indicate the fundamental and the first hve higher harmonics. 
Their spatial shapes agree well with the SPOD modes, which is not surprising as these 
modes have narrow spectral bands. Note however that each peak is split into several 
DFT modes, which indicates slight frequency variation. This becomes crucial for the 
higher harmonics, where the frequency jitter is significant and the mode energy is low. 
For the seventh harmonic, the DFT fails to reproduce the structure seen in the SPOD 
(figure 12(b)). The frequency variations detected by the SPOD are simply to high and 
the mono-frequent energy content too low. This emphasizes the superior noise rejection 
of the SPOD. 

In this example, the SPOD is superior to the POD and the DFT. The energy-ranked 
POD modes primarily suffer from the incompleteness of the data set. This is of immense 
importance as the relevant domain size is rarely known prior to a set of experiments, or 
POD analysis. The frequency-sharp DFT is insensitive to the domain size, but it fails to 
reconstruct weak modes with substantial noise and frequency jitter. The soft frequency 
constraint of the SPOD filter operation combines the advantages of both methods and 
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Figure 12. Fluidic oscillator: Results from SPOD for different filter lengths (a) Nf = 0 (POD), 
(b) Nf = 30 (SPOD), and (c) Nf — 2000 (DFT). For every filter length the SPOD spectrum 
is displayed as scatter plot (left), where a single dot indicates one mode pair (size and color 
Cij in (2.32)). For three selected pairs the spatial modes (upper row) and PSD of the temporal 
coefficient (lower row) are depicted. They are indicated by numbers in the SPOD spectrum, as 
well as between the small mode plots. The centerline of is indicated by a dash dotted line. 





























































































































22 


M. Sieher et al 



-2 0 2 -2 0 2 -2 0 2 
o-i/y/^ ^ily/1^ 

Figure 13. Phase portraits of first temporal coefficients from (a) POD, (b) SPOD, and (c) of 

both methods against each other. 

generates a clear mode space. The SPOD generates modes with distinct frequencies and 
mode shapes for modes even weaker than the overall noise level. 

For the fluidic oscillator, the DFT modes are nearly as accurate as those derived 
with the SPOD. The advantages of the SPOD are more obvious for less mono-frequent 
flow dynamics (see the previous examples). However, an additional advantage of the 
SPOD against the DFT is that it provides a reliable estimate of the oscillatory phase by 
accounting for the frequency jitter. Similar approaches, which also produce satisfactory 
results for the current case are described by Ostermann et al. (2015o), but again, the 
scope of SPOD is beyond this particular application. Figure 13(a,b) shows the phase 
portraits (Lissajous figures) of the temporal coefficients of the two most energetic POD 
(oi) and SPOD (6^) modes, respectively. The trajectory of the POD coefficients does 
not follow a clear circle that would indicate the limit-cycle. It rather follows a third of 
a circle and then collapses at one point. The coefficient of the SPOD modes follows a 
clear circle and the instantaneous phase and instantaneous frequency can consistently be 
deduced. A comparison of the first mode coefficient from both methods is shown in figure 
13(c). It reveals that half of the period is cut out for the POD (oi). This corresponds to 
the sweeping jet leaving the measurement domain, where the energy-based POD “sees” 
no jet. The SPOD properly recovers the temporal dynamics over the entire oscillation 
period. Furthermore, note that the SPOD produces coefficients with smooth temporal 
dynamics, while the POD coefficients show rather erratic movements. This is particularly 
important for reduced order modeling (Luchtenburg et al. 2009), phase averaging, and 
extended POD (Boree 2003). Practically, most of the further processing is eased if there 
is less noise. 


4. Summary and conclusion 

4.1. Properties and capabilities of the proposed method 
The SPOD is introduced as an extension of the POD for time-resolved data. This novel 
method involves a filter operation on the diagonal elements of the snapshot correlation 
matrix. The procedure is closely related to the classic snapshot POD with a negligible 
increase of algorithmic complexity and numerical costs. The SPOD filter allows for a 
continuous fading from the energetic optimality of POD to the spectral purity of DFT. 
It is conceptualized in a general form, with the POD and the DFT as the limiting cases. 
The concept of SPOD was developed trough our experience with experimental data 
processing, and not from a constraint optimization problem. It arose from the desperate 
need for a method that applies to a wide range of turbulent flows at minimum user input. 
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Figure 14. Schematic describing the main properties of the SPOD for increasing filter strength 
(from left to right). The top row shows pseudo-color plots of the filtered correlation matrix 
matrix (S). The phase portraits of the corresponding first two modes ( 6 i and 62 ) that describe 
the dominant oscillations are shown below. The axes of the plots shown here are the same as for 
the plots in figure 1 and 13, respectively. The graphs are based on the data already presented 
in section 2.2 and the SPOD is calculated from 200 snapshots. 


The SPOD is motivated based on theoretical considerations, where it is interpreted as a 
short time linearization of the flow dynamics. 

The key feature of the SPOD is the smoothing of the diagonal elements of the correla¬ 
tion matrix. This filter operation is shown to constrain the growth rates, amplitudes and 
frequencies of the SPOD modes. By setting the filter width, one gains control over the 
spectral bandwidth of the single modes. When the filter is set to the maximum length, 
the modes are assumed to be strictly periodic and the SPOD converges to the DPT. 

The principle of the SPOD is graphically summarized in figure 14. The images in the 
first row show the hltered correlation matrix at different filter widths Nf. The images 
below depict the phase portraits of the leading two modes (compare figure 13). It is 
apparent that the increased diagonal similarity of the correlation matrix, that goes in 
hand with the increased hlter width, successively limits the temporal variations of the 
mode amplitude and frequency until a stable limit cycle is reached. 

The application of SPOD to the flow held of a swirl-stabilized combustor, an airfoil 
with Gurney hap and a huidic oscillator revealed different advantageous features of the 
SPOD in comparison to POD or DPT. Por every single case there exist other suitable 
methods, which may perform equally well as the SPOD, but none of them is as versatile 
as the proposed method. 

The main advantages can be summarized as follows: 

Separation of structures: The soft spectral constraint of the SPOD allows for a 
much better separation of individual huid dynamic phenomena into single modes, whereas 
POD or DPT mix or spread them among several modes. 
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Noise rejection: SPOD it is insensitive to noise and even recovers dynamics that are 
weaker than the overall noise level. 

Data completion: SPOD can eliminate degradation of temporal dynamics of par¬ 
tially recorded phenomena. 

Plain dynamics: The mode coefficients are smooth in time and they feature ad¬ 
justable variations of frequency and amplitude (set by the filter size). 

The characteristics of the SPOD modes ease further processing such as the identifica¬ 
tion of linked modes, comparisons against other simultaneously acquired measurements 
(phase averaging or extended POD) and the identification of reduced order models. The 
SPOD may also provide a better basis for modal representation of snapshots as input for 
a DMD, as pursued in section 2.5. 


4.2. Concluding remarks 

The SPOD has proven to be a reliable method to extract distinct phenomena from 
turbulent flows. It was not derived from purely mathematical considerations, but rather 
evolved from practical data processing. Nonetheless, the method has well defined upper 
and lower bounds and generates modes that can be easily interpreted. As shown in the 
considered examples, SPOD is a simple way to extract coherent structures from turbulent 
flows, where the POD or the DPT failed to provide valuable results. The SPOD constrains 
the spectral content, but retains the modal sparsity of the POD. 

There are certainly plenty of other cases where this new method will ease the iden¬ 
tification of hidden coherent structures. Its true benefit lies in the fact that only one 
assumption is made about the investigated flow dynamics, which is the filter timescale. 
This can also be understood as an inertia imposed on the mode dynamics, limiting the 
rate of change of the frequency and amplitude. The choice of this timescale can be assessed 
from the flow’s dominant frequency or convective timescale, as shown in this article. The 
authors hope that the SPOD will give access to new fluid dynamical phenomena and 
enriches the available methods. 
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Appendix A. The spatial correlation version of SPOD 

The original POD can either be calculated from a spatial or temporal correlation, 
which allows a computationally efficient calculation by restricting the size of the problem 
to the number of the snapshots or the number of grid points. Similarly, the SPOD has a 
spatial correlation counterpart, which is computationally more efficient if the number of 
snapshots is much larger than the number of grid points. This approach is slightly more 
complex and less intuitive than the snapshot version. Nevertheless, it is very valuable 
if long time series of few sensors are supposed to be decomposed into proper modes 
to perform an extended POD or to derive the phase of an oscillatory mode from the 
measurements. Assume a simultaneous multi point pressure measurement that shall be 
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decomposed with SPOD. This series is decomposed as 

N 

P{x^,t) =p{xi) +'^bs{t)'I's{xi), (Al) 

where the number of measured positions M is much smaller than the number of samples 
N. The number of samples may easily reach a million or more, which complicates the 
solution in terms of memory requirements for the composition of the temporal correlation 
matrix and in terms of computational time for the solution of the eigenvalue problem. 
Therefore the temporal correlation described in section 2.3 is not feasible in this case. 
Instead, the spatial correlation should be employed, as outlined in this section. The multi 
time shift correlation tensor for the spatial SPOD reads 

Si,j,k,i = JPixi,t - kAt)p(xj,t - lAt)dt (A2) 

= k,l = -Nf...Nf, 

where p = P — p is the fluctuating part of the pressure and gk are the filter coeffi¬ 
cients. For numerical implementation this is reshaped to a matrix such that = 

S{i+k*M),{j+i*M): but for the theoretical description the tensor notation is retained. The 
correlation tensor is decomposed in eigenvalues and eigenvectors, such that 

Nf M 

E E Si,j,k,i^sixj,Ti) = ias'I'sixi,Tk) ; Ml ^ ^ ^ MAf > 0., (A3) 

l=-Nf j=l 

where = kAt. The eigenvector iFg constitutes a discrete convolution filter, which is 
applied to the time series to obtain the mode coefficients 

Nf M I - 

bs{t) = ^ '^J^^s{Xi,Tk)p{Xi,t-Tk). (A 4) 

k=-Nf i=l ' 

The spatial mode is the central part (jk = 0) of the convolution filter Ps{xi) = 'Ps{xi,0)- 
The entire eigenvectors Ps can be understood as a data driven filter bank, which allows 
for decomposition of time series into modal contributions. It might be applied to a single 
sensor, where each mode represents a certain spectral band of the signal. The principal 
approach is comparable to the empirical mode decomposition (Huang et al. 1998), but 
the SPOD can also handle multiple sensors. In case of multiple sensors, it gives excellent 
result when the phase of a dominant oscillation has to be reconstructed from pressure 
measurements. The approach outlined in this section is similar to the multi time delay 
POD phase estimation pursued by Hosseini et al. (2015). 

In contrast to the snapshot version, the computational cost of the spatial version of 
SPOD scales with the filter size. It is only more efficient than the snapshot approach if 
M{2Nf + 1) <N. 


Appendix B. Properties of the SPOD modes 

In section 2.3 it was shortly mentioned that the spatial SPOD modes are no longer 
orthogonal, which is only part of the truth. If the spatial mode P (A 3) together with all 
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of the temporally shifted instances is considered, they are orthonormal 

N; M 

'll ^j{xk,Tl) = 6ij. (Bl) 

l=-Nf k=l 

The snapshot based calculation introduced in section 2.3, however, is restricted to the 
zero time delay (t = 0) part of the spatial mode. Furthermore, the decomposition of 
the data into modal contributions is only feasible for time independent spatial modes. 
This limitation to one part of the the spatial mode 'Psixi) = ^sixi,T = 0) introduces 
some imperfections. The selected modes for the decomposition are neither normal nor 
orthogonal 

1 " 

'1'j{xk) 7^ bij. (B 2) 

k=l 

The loss of normality is a fact, whereas the norm of the spatial modes gives further 
insights to the data set. The norm 


0 




1 

M 


M 




(B3) 


indicates how well a single mode is represented by the investigated data set. 

With the application of the filter (2.9), an idealized correlation matrix is constructed 
that delivers modes, which are more or less captured by the initial data set. This fact 
is reflected by the deviation of the mode norm Q from one. Consider for example the 
measurement of the sweeping jet. There, the fundamental mode is only partially captured 
as shown in figure 13. With the SPOD, the missing data are completed and a SPOD 
mode pair with equal energy levels jii is obtained. However, for the construction of 
the spatial modes the coefficients are projected onto the original data (2.12). There, 
the imperfect representation of one of the two modes re-enters the processing. For the 
sweeping jet’s leading mode pair the norm of one mode is clearly below the other, but 
they approximately add up to one. Therefore, the eigenvalues describe the idealized 
energy content of the single modes and the norm of the spatial mode Ci corrects the 
deficits in comparison to the actual data set. The limiting SPOD cases (POD and DFT) 
do not show this deficits. The POD modes are already normalized Ci = 1 and for the 
DFT, the modes pair perfectly, while the norm of these pairs (f, j) exactly add up to one 

(0 + Ci = !)• 


REFERENCES 

Adrian, R. J. & Moin, P. 1988 Stochastic estimation of organized turbulent structure: 
homogeneous shear flow. Journal of Fluid Mechanics 190 , 531-559. 

Bach, A. B., Berg, R., Pechlivanoglou, G., Nayeri, C. & Paschereit, C. O. 2015a 
Experimental Investigation of the Aerodynamic Lift Response of an Active Finite Gurney 
Flap. In 53rd AIAA Aerospace Sciences Meeting, AIAA SciTech . American Institute of 
Aeronautics and Astronautics. 

Bach, A. B., Lennie, M., Pechlivanoglou, G., Nayeri, C. N. & Paschereit, C. O. 2014 
Finite micro-tab system for load control on a wind turbine. Journal of Physics: Conference 
Series 524 , 012082. 

Bach, A. B., Pechlivanoglou, G., Nayeri, C. & Pasghereit, C. O. 20156 Wake Vortex 
Field of an Airfoil Equipped with an Active Finite Gurney Flap. In 53rd AIAA Aerospace 
Sciences Meeting, AIAA SciTech . American Institute of Aeronautics and Astronautics. 



Spectral proper orthogonal decomposition 


27 


Berkooz, G., Holmes, P. & Lumley, J. L. 1993 The Proper Orthogonal Decomposition in 
the Analysis of Turbulent Flows. Annual Review of Fluid Mechanics 25 , 539-575. 

Boree, J. 2003 Extended proper orthogonal decomposition: a tool to analyse correlated events 
in turbulent flows. Experiments in Fluids 35 (2), 188-192. 

Bourgeois, J. A., Noack, B. R. & Martinuzzi, R. J. 2013 Generalized phase average with 
applications to sensor-based flow estimation of the wall-mounted square cylinder wake. 
Journal of Fluid Mechanies 736 , 316-350. 

Durgesh, V. & Naughton, J. W. 2010 Multi-time-delay LSE-POD complementary approach 
applied to unsteady high-Reynolds-number near wake flow. Experiments in Fluids 49 (3), 
571-583. 

Gray, R. M. 2005 Toeplitz and Circulant Matrices: A review. Foundations and Trends in 
Communications and Information Theory 2 (3), 155-239. 

Holmes, P., Lumley, J. L., Berkooz, G. & Rowley, G. W. 2012 Turbulence, Coherent 
Struetures, Dynamical Systems and Symmetry, 2nd edn. Gambridge and UK and New 
York: Cambridge University Press. 

Hosseini, Z., Martinuzzi, R. J. & Noack, B. R. 2015 Sensor-based estimation of the velocity 
in the wake of a low-aspect-ratio pyramid. Experiments in Fluids 56 (1). 

Huang, N. E., Shen, Z., Long, S. R., Wu, M. C., Shih, H.H., Zheng, Q., Yen, N.-C., 
Tung, C. C. & Liu, H. H. 1998 The empirical mode decomposition and the hilbert 
spectrum for nonlinear and non-stationary time series analysis. Proceedings of the Royal 
Society of London A: Mathematical, Physical and Engineering Sciences 454 (1971), 903- 
995. 

Lasagna, D., Orazi, M. & luso, G. 2013 Multi-time delay, multi-point linear stochastic 
estimation of a cavity shear layer velocity from wall-pressure measurements. Physics of 
Fluids 25 (1), 017101. 

Luchtenburg, D. M., Gunther, B., Noack, B. R., King, R. & Tadmor, G. 2009 A 
generalized mean-field model of the natural and high-frequency actuated flow around 
a high-lift configuration. Journal of Fluid Mechanics 623 , 283-316. 

Lumley, J. L. 1970 Stochastic tools in turbulence, Applied mathematics and mechanics, vol. v. 
12. New York: Academic Press. 

Oberleithner, K., Rukes, L. & Soria, J. 2014 Mean flow stability analysis of oscillating jet 
experiments. Journal of Fluid Mechanics 757, 1-32. 

Oberleithner, K., Sieber, M., Nayeri, C. N., Paschereit, C. O., Petz, C., Hege, H.-C., 
Noack, B. R. & Wygnanski, I. 2011 Three-dimensional coherent structures in a swirling 
jet undergoing vortex breakdown: stability analysis and empirical mode construction. 
Journal of Fluid Mechanics 679 , 383-414. 

Oberleithner, K., Stohr, M., Seong, H. L, Arndt, C. M. & Steinberg, A. M. 2015 
Formation and flame-induced supression of the processing vortex core in a swirl combustor: 
Experiments and linear stability analysis. Combustion and Flame 162 (8), 3100-3114. 

OSTERMANN, F., WOSZIDLO, R., GAERTLEIN, S., NAYERI, G. & PASCHEREIT, C. O. 2015a 
Phase-Averaging Methods for the Natural Flowfield of a Fluidic Oscillator. AIAA Journal 
53 (8), 2359-2368. 

OSTERMANN, F., WosziDLO, R., Nayeri, C. & PASCHEREIT, C. O. 2015& Experimental 
Comparison between the Flow Field of Two Common Fluidic Oscillator Designs. In 53rd 
AIAA Aerospace Sciences Meeting, AIAA SciTech . American Institute of Aeronautics 
and Astronautics. 

Raffel, M., Kompenhans, J., Wereley, S. T. & Willert, C. E. 2007 Particle image 
velocimetry: A practical guide, 2nd edn. Heidelberg and New York: Springer. 

Raiola, M., Discetti, S. & Ianiro, A. 2015 On PIV random error minimization with optimal 
POD-based low-order reconstruction. Experiments in Fluids 56 (4). 

Reichel, T. G., Terhaar, S. & Paschereit, O. 2015 Increasing Flashback Resistance in 
Lean Premixed Swirl-Stabilized Hydrogen Combustion by Axial Air Injection. Journal of 
Engineering for Cas Turbines and Power 137 (7), 071503. 

Rowley, C. W., Mezic, L, Bagheri, S., Schlatter, P. & Henningson, D. S. 2009 Spectral 
analysis of nonlinear flows. Journal of Fluid Mechanics 641 , 115-127. 

SCARANO, F. 2002 Iterative image deformation methods in PIV. Measurement Science and 
Teehnology 13 (1), 1-19. 



28 


M. Sieber et al 


SCHLEGEL, M., NOACK, B. R., JORDAN, P., DlLLMANN, A., GROSCHEL, E., SCHRODER, W., 
Wei, M., Freund, J. B., Lehmann, O. & Tadmor, G. 2012 On least-order flow 
representations for aerodynamics and aeroacoustics. Journal of Fluid Mechanics 697 , 
367-398. 

Schmid, P. J. 2010 Dynamic mode decomposition of numerical and experimental data. Journal 
of Fluid Mechanics 656 , 5-28. 

SiROViCH, L. 1987 Turbulence and the dynamics of coherent structures: Part I: Coherent 
Structures. Quarterly of Applied Mathematics 45 (3), 561-571. 

Terhaar, S. 2015 Identification and modeling of coherent structures in swirl stabilized 
combustors at dry and steam diluted conditions. PhD thesis, Technische Universitat 
Berlin. 

Terhaar, S., Oberleithner, K. & Paschereit, C. O. 2015 Key parameters governing the 
precessing vortex core in reacting flows: An experimental and analytical study. Proceedings 
of the Combustion Institute 35 (3), 3347-3354. 

Troolin, D. R., Longmire, E. K. & Lai, W. T. 2006 Time resolved PIV analysis of flow over 
a NACA 0015 airfoil with Gurney flap. Experiments in Fluids 41 (2), 241-254. 

Tu, J. H., Rowley, C. W., Luchtenburg, D. M., Brunton, S. L. & Kutz, J. N. 2014 
On dynamic mode decomposition: Theory and applications. Journal of Computational 
Dynamics 1 (2), 391-421. 

WiLLERT, C. E. & Gharib, M. 1991 Digital particle image velocimetry. Experiments in Fluids 
10 (4). 

Williams, M. O., Kevrekidis, I. G. & Rowley, C. W. Submitted A data-driven 
approximation of the Koopman operator: extending dynamic mode decomposition. arXiv 
preprint arXiv:1408.4408 . 

Wise, J. 1955 The autocorrelation function and the spectral density function. Biometrika 42 (1- 
2), 151-159. 

Woszidlo, R., Ostermann, F., Nayeri, C. N. & Pasghereit, C. O. 2015 The time-resolved 
natural flow field of a fluidic oscillator. Experiments in Eluids 56 (6). 

Woszidlo, R., Stumper, T., Nayeri, G. & Pasghereit, C. O. 2014 Experimental study 
on bluff body drag reduction with fluidic oscillators. In 52rd AIAA Aerospace Sciences 
Meeting, AIAA SciTech . American Institute of Aeronautics and Astronautics. 



